Open In Colab | Descargar notebook.

Instalar librerias necesarias (en modo silencioso -q):

[1]:

!wget -q https://raw.githubusercontent.com/arqlm/arqlm.github.io/main/_static/libraries.txt
!pip install -q -r /content/libraries.txt
!pip install -q --upgrade plotly[kaleido]
!rm -r /content/libraries.txt
!rm -r /content/sample_data/ ## Esta linea no es necesaria cuando no se trabaja en colab.google

   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.9/7.9 MB 35.8 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 51.5/51.5 kB 1.7 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 9.6/9.6 MB 64.7 MB/s eta 0:00:00
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 51.3/51.3 kB 2.9 MB/s eta 0:00:00


[29]:

## Importar librerias

import pandas as pd
import numpy as np
import base64
import datetime
import io

import plotly.io as pio
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode

1. Introducción

Tan pronto como comenzamos a observar un conjunto de muestras para inferir propiedades que les son comunes, entramos en el ámbito de la inferencia. Nos interesa encontrar el conjunto de muestras y el volumen asociado, donde las variables se comportan de manera consistente. Hasta ahora, no hemos definido claramente esta consistencia, pero digamos por ahora que queremos que las variables mantengan las mismas propiedades dentro de este volumen. Esto requiere una abstracción para ir más allá de los datos. Sería imposible aprender sobre una ubicación no muestreada, a menos que asumamos algo sobre las propiedades de las variables de interés en esa ubicación. El paso natural es asumir que “se comporta de manera similar a sus vecinos”. Cuando decimos esto, nos referimos a dos aspectos diferentes: primero, que las propiedades estadísticas son similares, y segundo, que las propiedades espaciales también son similares a las de los vecinos.

1.1. Una noción sobre lo que definiremos como Función Aleatoria

El valor desconocido en una ubicación no muestreada puede modelarse utilizando el concepto de variable aleatoria. Esto significa que, en cada ubicación, la variable regionalizada \(z(\mathbf{u})\) se modela como una variable aleatoria \(Z(\mathbf{u})\) que toma valores de acuerdo a una distribución de probabilidad. Nótese la notación: minúscula para la variable medible real y mayúscula para el objeto matemático. El objetivo principal de los pasos que se presentarán en las siguientes secciones y capítulos es inferir las características de la variable aleatoria en ubicaciones no muestreadas.

Además de este objetivo perfectamente razonable de inferencia local, también nos interesará comprender la relación entre las diferentes ubicaciones. Las variables aleatorias \(Z(\mathbf{u}_1)\) y \(Z(\mathbf{u}_2)\) están relacionadas porque existe correlación espacial. Podemos extender esta noción a la relación entre múltiples puntos. El conjunto de variables aleatorias dentro de un dominio se conoce como función aleatoria. La función aleatoria describe las propiedades estadísticas de las variables aleatorias dentro del dominio y sus relaciones espaciales. En geoestadística, la inferencia está orientada a caracterizar el comportamiento de la función aleatoria. Veremos que esto requiere un análisis tanto estadístico como espacial.

Definiciones:

  • Variable regionalizada: es la variable que puede medirse en el mundo real.

  • Variable aleatoria: es una abstracción matemática, donde se asume que el valor desconocido en una ubicación no muestreada sigue una distribución de probabilidad.

  • Función aleatoria: es el conjunto de variables aleatorias en un dominio definido y debe caracterizarse desde un punto de vista estadístico y espacial.

  • Dominio: es el volumen dentro del cual las variables aleatorias tienen un comportamiento homogéneo (estacionario). Nótese que la estacionariedad es una decisión tomada con fines de modelado.

Procedamos a continuación a la búsqueda de controles geológicos que nos ayuden a definir mas tarde dominios de estimación.

2. Búsqueda de controles geológicos (por Zona Mineral)

Continuando con el ejemplo presentado anteriormente, analizaremos las distribuciones para determinar dominios razonables para la estimación y simulación.

2.1. Carga de archivos

[30]:

import pandas as pd
import numpy as np

import warnings
warnings.simplefilter(action="ignore")

# Cargar base de datos en pandas
DH = pd.read_csv('EvYacData.csv', sep=',', encoding='latin1')
DH['cu_pct'][DH['cu_pct'] <= 0] = np.nan

DH

[30]:
Este Norte Elevación au_ppm ag_ppm cu_pct aucn_ppm cucn_ppm Zmin Alte Lito
0 472186.686 6925804.447 4220.763 -99.00 -99.0 NaN -99.0 -99.0 OXI ILL_CLO IBX_MM
1 472187.202 6925805.493 4213.861 0.30 2.3 0.011 -99.0 -99.0 OXI ILL_CLO IBX_MM
2 472187.343 6925805.770 4211.986 0.47 16.2 0.032 -99.0 -99.0 OXI ILL_CLO IBX_MM
3 472187.493 6925806.060 4210.013 0.31 2.3 0.018 -99.0 -99.0 OXI ILL_CLO IBX_MM
4 472187.642 6925806.351 4208.040 0.29 2.1 0.010 -99.0 -99.0 OXI ILL_CLO IBX_MM
... ... ... ... ... ... ... ... ... ... ... ...
77836 471293.198 6925823.101 3761.674 0.03 1.0 0.002 -99.0 -99.0 BACK PROP ECS
77837 471293.540 6925824.040 3759.942 0.01 0.7 0.001 -99.0 -99.0 BACK PROP ECS
77838 471293.882 6925824.980 3758.209 0.01 0.4 0.001 -99.0 -99.0 BACK PROP ECS
77839 471294.224 6925825.920 3756.477 0.06 1.3 0.004 -99.0 -99.0 BACK PROP ECS
77840 471294.607 6925826.972 3754.538 0.06 0.5 0.003 -99.0 -99.0 BACK PROP ECS

77841 rows × 11 columns

Si observamos la distribución espacial de las muestras, codificadas según zonas mineral, podemos ver la ubicación espacial de las diferentes zonas minerales (Figura 1).

[31]:

DH["Zmin"] = DH["Zmin"].astype(str)
fig = px.scatter_3d(DH, x='Este', y='Norte', z='Elevación', color='Zmin',category_orders={'Zmin': sorted(DH.groupby('Zmin').groups.keys())})
fig.update_traces(marker=dict(size=2.0))

# Mejor estilo para leyenda
fig.update_layout(
    hovermode=False,
    legend_title_text='Zona Mineral',
    legend=dict(
        itemsizing='constant',
        font=dict(size=14),
        bgcolor='rgba(255,255,255,0.8)',
        bordercolor='black',
        borderwidth=1,
    )
)

fig.show()

El ejemplo presentado anteriormente corresponde a un caso típico de mineralización estratificada, donde las zonas minerales (por ejemplo, OXIdes - Zona Oxidada, TRANSition, BACKground - Zona Primaria) se distribuyen en capas o estratos bien definidos a lo largo del perfil geológico. Este tipo de mineralización es común en depósitos donde los procesos geoquímicos y geológicos han generado una zonificación vertical.

En estos sistemas, la mineralización se organiza en bandas o capas superpuestas, reflejando la alteración progresiva de los minerales primarios y la migración de elementos a lo largo del tiempo. La visualización espacial de las muestras, codificadas por zona mineral, permite identificar claramente estos dominios y es fundamental para la definición de dominios de estimación y simulación geoestadística.

A continuación se muestra un esquema típico de un perfil de mineralización estratificada, tomado como referencia de archaeometallurgie.de:

Perfil típico de mineralización estratificada

Referencia:
https://en.archaeometallurgie.de/gossan-iron-cap/

2.2. Estacionariedad

La estacionariedad se refiere a una noción de homogeneidad dentro de un dominio. Intuitivamente, esperamos encontrar el mismo comportamiento en cualquier parte de un dominio. Cuando hablamos de comportamiento, nos referimos a propiedades estadísticas y propiedades espaciales.

Por ejemplo, un dominio homogéneo debería tener el mismo valor medio para la variable en diferentes partes de ese dominio. No esperaríamos encontrar un valor homogéneo, es decir, el mismo valor en cada lugar del dominio, pero sí esperamos encontrar valores dentro del mismo rango. De manera similar, si los valores muestran un cierto nivel de variabilidad en una parte del dominio, esperamos encontrar la misma variabilidad en otras partes del dominio.

La estacionariedad puede clasificarse como:

  • De primer orden.

  • De segundo orden.

  • Cuasi de segundo orden.

  • Estacionariedad estricta.

La estacionariedad de primer orden implica que las medias de las variables aleatorias dentro del dominio de trabajo \(D\) son constantes: \begin{equation} m(Z(\mathbf{u}))= costante \qquad \forall \mathbf{u} \in D \end{equation}

Esto podría verificarse a partir de los datos disponibles calculando un promedio ya seas por categoria geologica (tipo de roca o litología, alteración, zona mineral), o bien dentro de una ventana móvil en caso de no tener disponibles estas categorias. Si se observa un cambio sistemático en el valor promedio a lo largo del dominio, podríamos decir que no se cumple la estacionariedad de primer orden.

La estacionariedad de segundo orden requiere, además de la condición sobre las medias, que los momentos de segundo orden de las variables aleatorias sean constantes. Esto significa que las varianzas locales son relativamente constantes en el dominio:

\begin{equation} \sigma^2(Z(\mathbf{u}))= costante \qquad \forall \mathbf{u} \in D \end{equation}

También significa que las relaciones de dos puntos son constantes en todo el dominio. Así, deberíamos esperar la misma cantidad de similitud entre puntos separados por una distancia \(\mathbf{h}\) (que en realidad es un vector, con magnitud y dirección) en cualquier posición dentro del dominio \(D\). Volveremos a esto más adelante, una vez que introduzcamos las medidas de continuidad espacial.

La estacionariedad cuasi de segundo orden es una aproximación conveniente del tipo de segundo orden. Básicamente, relaja las condiciónes anteriormente vistas a vecindarios locales \(\mathcal{V}(\mathbf{u})\) alrededor de \(\mathbf{u}\):

\begin{equation} m(Z(\mathbf{u}))= costante \qquad \forall \mathbf{u} \in \mathcal{V}(\mathbf{u}) \\ \sigma^2(Z(\mathbf{u}))= costante \qquad \forall \mathbf{u} \in \mathcal{V}(\mathbf{u}) \end{equation}

Por ejemplo, podemos decir que la función aleatoria es estacionaria de segundo orden solo dentro del vecindario \(\mathcal{V}(\mathbf{u})\) alrededor de una ubicación \((\mathbf{u})\) que estamos tratando de caracterizar. Esto realmente simplifica la decisión de estacionariedad, ya que puede existir un cambio en la media local (o en la variabilidad), pero aún podemos asumir que dentro de cada vecindario local las cosas permanecen constantes. Esta decisión se utiliza más adelante en el *kriging ordinario**.

La estacionariedad estricta es la suposición de que todos los momentos de la función aleatoria permanecen constantes en todo el dominio. Los momentos se refieren a diferentes estadísticas de la función aleatoria, como la media, la varianza y la continuidad espacial medida por estadísticas de dos puntos, o por estadísticas de patrones, es decir, estadísticas que relacionan múltiples puntos al mismo tiempo. Esta es una suposición poco realista, pero puede ser necesaria para la inferencia y para desarrollar algunas herramientas teóricas.

En resumen, cualquier cambio sistemático obvio en las propiedades de la función aleatoria a lo largo del dominio debe verse como una advertencia en cuanto a las suposiciones que necesitamos hacer sobre la función aleatoria cuando realizamos inferencias sobre diferentes estadísticas.

Procedamos a continuación a verificar que hipótesis podremos asumir basándonos en las estadísticas por categoría geológica.

2.3. Distribución estadística por Zona Mineral

Calculamos los gráficos de probabilidad de la ley de cobre para los diferentes tipos de zona mineral (Figura 2). La descripción de los tipos de roca se presenta en la Tabla 1.

[32]:

import numpy as np
import plotly.graph_objects as go
import scipy.stats as stats
from statsmodels.distributions.empirical_distribution import ECDF

# Plot all probability curves in the same figure
fig = go.Figure()
for category in sorted(DH.groupby('Zmin').groups.keys()):
    values = DH.groupby('Zmin').get_group(category)['cu_pct'].dropna().values
    # Sort values
    values_sorted = np.sort(values)
    n = len(values_sorted)
    cumprob = (np.arange(1, n+1) - 0.5) / n
    normal_scores = stats.norm.ppf(cumprob)
    fig.add_trace(go.Scatter(
        x=values_sorted,
        y=normal_scores,
        mode='markers',
        marker=dict(size=5, opacity=0.7),
        name=str(category)
    ))

# Custom y-axis to match imagen.png
y_ticks = stats.norm.ppf([0.0001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.99999])
y_ticklabels = ['0.001%','1%', '5%', '10%', '20%', '30%', '40%', '50%', '60%', '70%', '80%', '90%', '95%', '99%', '99.999%']

fig.update_layout(
    title="Gráfico de probabilidad log-normal por Zona Mineral",
    title_x=0.5,
    xaxis_title='Cu [%]',
    xaxis=dict(
        type='log',
    ),
    yaxis_title='Probabilidad',
    yaxis=dict(
        tickmode='array',
        tickvals=y_ticks,
        ticktext=y_ticklabels,
        range=[stats.norm.ppf(0.00001), stats.norm.ppf(0.99999)]
    ),
    width=800,
    height=600,
    font_family="Times New Roman",
    font_size=14,
    font_color="black",
    plot_bgcolor='white',
    autosize=False,
    showlegend=True
)
fig.update_layout(
    legend_title_text='Zona Mineral',
    legend=dict(
        itemsizing='constant',
        font=dict(size=14),
        bgcolor='rgba(255,255,255,0.8)',
        bordercolor='black',
        borderwidth=1
    )
)
fig.update_xaxes(gridcolor='rgba(0,0,0,0.1)')
fig.update_yaxes(gridcolor='rgba(0,0,0,0.1)')
fig.update_layout(hovermode=False)

fig.show()


[33]:

stats_by_mine = DH.groupby(['Zmin']).describe().round(2)
stats_by_mine['cu_pct']

[33]:
count mean std min 25% 50% 75% max
Zmin
0 3030.0 0.01 0.02 0.0 0.00 0.00 0.00 0.27
BACK 46632.0 0.19 0.22 0.0 0.05 0.15 0.27 8.30
OXI 17266.0 0.06 0.10 0.0 0.02 0.03 0.05 1.82
TRANS 9979.0 0.12 0.21 0.0 0.02 0.04 0.15 4.28

Es posible identificar que la Unidad 20 corresponde a la BACK (en Rojo), la unidad más relevante, con mayores leyes y volumen dentro del yacimiento. Las otras unidades: óxidos y zona de transición, le sigue en numero de muestras, presentando leyes más bajas y rodeando a la zona de background. La zona sin definir, 0, es una zona prácticamente sin leyes. Para verificar las distintas distribución de leyes, podemos desplegar las muestras por separado, de acuerdo a su código por zona mineral.

[34]:

DH_filtered = DH[(DH['Zmin'] == 'TRANS') | (DH['Zmin'] == 'BACK')] # Agrupamos las dos zonas con mayor ley
fig = px.scatter_3d(DH_filtered, x='Este', y='Norte', z='Elevación', color='cu_pct', color_continuous_scale=px.colors.sequential.Rainbow[1:], range_color=[0.0, DH['cu_pct'].quantile(0.95)])
fig.update_traces(marker=dict(size=2.0))
fig.update_layout(hovermode=False)
fig.show()

DH_filtered = DH[(DH['Zmin'] == 'OXI')]
fig = px.scatter_3d(DH_filtered, x='Este', y='Norte', z='Elevación', color='cu_pct', color_continuous_scale=px.colors.sequential.Rainbow[1:], range_color=[0.0, DH['cu_pct'].quantile(0.95)])
fig.update_traces(marker=dict(size=2.0))
fig.update_layout(hovermode=False)
fig.show()

Este ejercicio se puede repetir para las otras dos características geológicas restantes: Alteración y Litología